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Abstract. The IceCube Neutrino Observatory has provided the first map of the high energy 
(~ 0.01 - 1 PeV) sky in neutrinos. Since neutrinos propagate undeflected, their arrival direc¬ 
tion is an important identifier for sources of high energy particle acceleration. Reconstructed 
arrival directions are consistent with an extragalactic origin, with possibly a galactic compo¬ 
nent, of the neutrino flux. We present a statistical analysis of positional coincidences of the 
IceCube neutrinos with known astrophysical objects from several catalogs. When considering 
starburst galaxies with the highest flux in gamma-rays and infrared radiation, up to n = 8 
coincidences are found, representing an excess over the ~ 4 predicted for the randomized, or 
“null” distribution. The probability that this excess is realized in the null case, the p- value, is 
p = 0.042. This value falls to p = 0.003 for a partial subset of gamma-ray-detected starburst 
galaxies and superbubble regions in the galactic neighborhood. Therefore, it is possible that 
starburst galaxies, and the typically hundreds of superbubble regions within them, might 
account for a portion of IceCube neutrinos. The physical plausibility of such correlation is 
discussed briefly. 

Keywords: neutrino astronomy, neutrino experiments, star formation, gamma ray experi¬ 
ments 


ArXiv ePrint: 1507.05711 


Contents 


1 Introduction 1 

2 Data 3 

2.1 IceCube neutrino detections 3 

2.2 Catalogs of possible counterparts 3 

3 Statistical method 4 

3.1 Combining neutrino events 6 

4 Analysis of possible counterparts 6 

4.1 “Null” results 7 

4.2 Star-forming activity 7 

4.2.1 Star burst galaxies 8 

4.2.2 Starburst galaxies and local superbubble regions 11 

5 Discussion: star formation as a possible origin of high energy neutrinos 12 

6 Summary and conclusion 16 

A Formalism: the null hypothesis 18 

B Null results 19 

1 Introduction 


Extra-solar neutrino astronomy is an infant science, born in 2013, with the first detection of 
astrophysical neutrinos of energies up to ~PeV at the IceCube experiment in Antarctica [1-3] . 
The origin of these neutrinos has not been established yet and represents an important goal 
to learn about the fundamental physics at play in astrophysical accelerators. The inherent 
properties of neutrinos - neutral, weakly interacting - offer unique probes into relatively 
unknown high energy mechanisms such as stellar core collapse and jet formation, particle 
acceleration in magnetic fields, shockwave propagation, etc. 

Theoretically, there is a close relationship between neutrinos and cosmic ray protons 
(CRp). Comparable fluxes of neutrinos and gamma rays are expected as by-products of CRp 
interactions. Neutrinos can be created in proton-proton (pp) interactions and subsequent 
cascades of charged and neutral pions: 

p + p —> 7T° + + anything 

7 T° —> 7 T 7 
-» fJ^ + 

1^ ->• e ± + v e {v e ) + (1-1) 

or in proton-photon (jyy) interactions, e.g., with cosmic microwave background (CMB) radi¬ 
ation: 

p + 7 —t n + n + /p + 7T° (1.2) 
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with subsequent pion decays as in eq. (1.1). Neutrinos created in the above interactions have 
~ 5% of the initial CRp energy and ~ 50 — 75% of the gamma-ray energy [3] . 

In light of the neutrino-gamma-ray connection, the search for sources of the IceCube 
neutrinos has turned to the most powerful known gamma-ray emitters, and, more broadly, 
to the objects that show high energy activity. Since the IceCube data are consistent with 
a diffuse flux, several authors have compared their energy distribution with the spectra 
predicted theoretically for the diffuse neutrino flux for several classes of possible sources. 
It was found that the pp mechanism appears to naturally fit the data (see e.g., [4]), and, 
among the pp-dominated objects, starburst and star-forming galaxies have emerged as a 
particularly interesting possibility, fitting the data well both in spectrum and normalization 
[5-8]. In these objects, both the elements required for abundant neutrino production - 
proton acceleration and a dense proton background - are expected. Indeed, starburst and 
star-forming galaxies are characterized by their high rate of star formation, which implies a 
high rate of proton-accelerating jets in core collapse supernovae and/or supernova remnants. 
Because star formation typically occurs in dense hydrogen clouds, these galaxies should also 
be good proton absorbers, and thus efficient neutrino sources [9-11]. 

In parallel, searches have been conducted for positional associations of the neutrino 
data to specific objects. This task is made difficult by the poor angular resolution of Ice- 
Cube (~15° [2]). However, positional matching is attractive because it is practically model- 
independent, relying only on the fact that - in the absence of exotic physics - neutrinos 
propagate undeflected from the production point to Earth. The searches performed by the 
IceCube collaboration, including point-like and extended-emission sources, all gave results 
consistent with background only [3, 12, 13]. Other authors have pointed out non-significant 
associations of some of the data with galactic objects, mainly the Galactic Center [14, 15] and 
Fermi bubbles [14, 16]. Coincidences with up to 93% confidence level have also been noted 
with the arrival directions of ultra-high energy cosmic rays (UHECR) [17, 18]. Several classes 
of extragalactic point-sources have been examined as well for spatial associations, in partic¬ 
ular blazars [19-26], Seyfert galaxies [18] and star-forming galaxies with high luminosities in 
hydrogen cyanide (HCN) emission [27]. The conclusions were mixed in these cases. 

At this time, the status of searches for positional associations of the IceCube data with 
astrophysical objects is still heterogeneous, and no consistent picture has emerged. There 
remains a need for advancement towards a more systematic and interdisciplinary approach, 
that can fully incorporate knowledge and techniques from both neutrino physics and as¬ 
tronomy. This approach could be used to test the hypothesis that high energy neutrinos 
ultimately originate from star formation. 

This paper is meant to be a step in this direction. We perform a statistical analysis 
of the IceCube neutrino data, to test for spatial coincidences with the brightest candidate 
sources from several catalogs. Compared to previous literature, our work has a stronger focus 
on star formation as a possible origin of the neutrino events. It is on this subject that our 
results are the most significant. 

Section 2 proceeds with a description of the IceCube neutrino data set and the catalogs 
that we have considered for counterparts. In section 3, we discuss the statistical method used. 
The motivation and results of the statistical analyses focusing on star-forming activity are 
presented in section 4, followed by a discussion of the possible role of starburst galaxies and 
nearby star-forming regions in section 5. Additional testing of candidate sources for which 
we found null results has been included in appendix B. Conclusions are given in section 6. 
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Where relevant, we use the Planck 2015 cosmological parameters (Hq = 68 km s 1 Mpc , 
n M = 0.32, n A = 0.68) [28], 

2 Data 

2.1 IceCube neutrino detections 

We consider the 37 data events obtained by IceCube after 988 days of running [3]. When 
needed, individual events will be referred to by their number as in Table 1 of ref. [3]. Each 
event is classified as track-like or shower-like, depending on its topology in the detector. 
The track-like events occur when a neutrino interaction results in a particle shower with 
a discernible muon track and therefore have smaller angular resolution on the sky (< 1°). 
The nine observed track-like events are consistent with the expected background of 8.4 ± 4.2 
atmospheric muons. Shower-like events, on the other hand, result in a spherical light pattern 
produced by particle showers with no discernible muon and therefore have poorer angular 
resolution (median ~15°; 50% confidence level of positional errors). The twenty-eight shower¬ 
like events are in excess of the expected background of 6.6^(] atmospheric neutrinos. Events 
28 and 32 have coincident detections at the Ice Top surface array. Thus, they have been 
identified as background [3] and will be excluded from our analysis. 

The 35 events we use are shown in e.g., figure 1, in equatorial J2000 coordinates, with 
their median angular errors. Because IceCube is located near the South Pole, its horizon 
coincides with the celestial Equator. Due to absorption of neutrinos in the Earth, the detector 
is considerably less sensitive to up-going neutrinos (i.e. below the horizon, as coming from 
the northern sky) compared with down-going neutrinos (i.e. above the horizon, from the 
southern sky). The difference in performance increases with increasing neutrino energy — at 
E ~ PeV, the Earth is essentially opaque. This feature explains the noticeable asymmetry 
in the event distribution between the two hemispheres as seen in e.g., figure 1. 

2.2 Catalogs of possible counterparts 

When searching astronomical catalogs for possible counterparts of the IceCube neutrinos, 
it is logical to consider gamma-ray emitters detected in the same energy range (E > 100 
TeV) as the IceCube events. However, gamma-rays with E 7 > 100 GeV can suffer from 
absorption due to photon interactions with extragalactic background light and at E 7 > 
PeV, interactions with the cosmic microwave background [29]. Furthermore, current TeV 
observations lack uniform and complete all-sky coverage, which is one of the conditions of 
validity of our analysis (see section 3). Therefore, we resort to primarily using observations 
from the Fermi Large Area Telescope (hereafter, Fermi- LAT) [30] with sensitivity up to 
~ 300 GeV, and specifically the Fermi -LAT 3FGL catalog [31]. We then supplement the 
analysis with TeVCat detections and a catalog of starburst galaxies based on their infra-red 
(IR) flux [32, 33], 

Considering that many modern, all-sky catalogs include thousands of objects, it was 
necessary to apply selection criteria to each catalog. Two main principles are used to choose 
sensible selection criteria. The first is uniformity: each set of candidates is made of sources 
of the same type/morphology (e.g., blazars, Seyfert galaxies, etc.). The second principle is 
the assumption that, in a given class of viable candidates, those that appear the brightest 
in photon flux should also be the brightest in neutrino flux, and therefore most likely to 
be responsible for the neutrino events. Hence, lower limits of photon flux at appropriate 
wavelengths will be imposed. 
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Finally, we emphasize that our selection procedure is completely blind with respect to 
the position of a neutrino source candidate in the sky. 

We investigated the following catalogs. 

TeVCat. TeVCat 1 is a compilation of currently and previously known TeV gamma-ray 
sources. TeV gamma-ray instruments, sensitive to energies between 100 GeV - 100 TeV, can 
image ganrma-ray emission via atmospheric Cherenkov telescope arrays. The procedure of 
reconstructing particle showers, created from the interaction of gamma-rays with the Earth’s 
atmosphere, allows an angular resolution of photon arrival typically of < 0.1° at 1 TeV 
[34, 35], 

3FGL. Fermi- LAT is a pair-production ganrma-ray instrument operating in the 20 MeV 
- 300 GeV range, and the current, main workhorse of space based gamma-ray observations. 
Its angular resolution varies from ~5° at 100 MeV to 0.8° at 1 GeV. The LAT 4-year Point 
Source Catalog (3FGL) covers the entire sky for at least 15 Ms of observing time, reaching 
a detection threshold of ~ 3 x 10” 12 erg cm” 2 3 s” 1 in the 100 MeV — 100 GeV energy range 
[31]. It is known that at low galactic latitudes (|b| < 10°) diffuse emission from pion decay, 
bremsstrahlung, and inverse Compton scattering reduces the sensitivity, although it is hard 
to quantify by exactly how much [31]. 

IRAS. We use flux measurements in the mid to far IR as an indicator of star-forming 
activity in starburst galaxies [32], In particular, we chose objects with the highest fluxes at 
100 fi m, which is close to the spectral luminosity peak of heated dust in starburst galaxies. 
The 100 fi m measurements were all gathered with the Infrared Astronomical Satellite (IRAS) 
Revised Bright Galaxy Sample [33], a sample of galaxies chosen to have 60 fim flux density 
greater than 5.24 Jy, 2 which covers 93% of the sky excluding only a strip within 5° of the 
galactic plane (|6| < 5°). 

3 Statistical method 

We adopt a version of the likelihood ratio statistical method which is commonly used in 
astronomy (see, e.g. [36-38]) and has been used in high energy astrophysics [18, 39] to 
test the spatial correlation between a set of data points (the neutrino IceCube data) and a 
population of candidate sources. The statistical variable of interest is the angular distance 
between each neutrino event and the candidate source closest to it. If the data and the 
potential sources are causally related, we expect an abundance in low distances. This forms 
the basic premise of our analysis. 

First, we define a dimensionless distance which is weighted by statistical spatial errors. 
Consider two objects in the sky with equatorial coordinates (oti,Si) and ( aj,Sj ), with the 
first coordinate being the Right Ascension (RA) and the second the declination (dec), and 
angular positional errors cq and (Jj? Their angular separation is then: 

Sij = cos” 1 (sin (Si) sin (Sj) + cos (Si) cos (Sj) cos (Acc^)), (3-1) 

1 http://tevcat.uchicago.edu/ 

2 1 Jy (Jansky)= 10” 23 erg s” 1 cm” 2 Hz” 1 . 

3 We assume that positional errors are symmetric in all directions, resulting in a spherical cone around each 
measured position. This is the case for the neutrino events considered here. 
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where Actjj is the difference in right ascension coordinates, and their weighted, dimensionless 
angular distance can be defined as: 


Rij — 


S\ 


*3 


a i + a j 


(3.2) 


The next step is to consider the set of i = 1,2,. ,N (N = 35) neutrino data and j = 

1,2 M candidates of a certain class, and, for each datum i. find the distance to the 
closest candidate 4 : 

n = Min (jjRij . (3.3) 

We then have N values of r - the index i will be dropped from here on out for simplicity 
of notation. In all cases considered in this work, the angular errors of the IceCube data 
dominate over the positional uncertainty of the sources, which therefore are neglected i.e., 
(Tj =0. If a neutrino event has r < 1, its positional error encompasses the nearest candidate’s 
on-sky location. In what follows, this condition will be used as an indicator of a plausible 
positional correlation between the datum and the candidate. Of particular interest will be 
the number of data for which the weighted distance is r < 1. 

The final steps consist of generating the distribution of the variable r and comparing 
it with a null distribution. The latter is obtained through the hypothesis of a uniform 
distribution of sources in the sky. It represents the outcome expected if the data and the 
candidate sources are not causally related, and spatial coincidences are simply the result 
of random accident. The null distribution can be calculated, using a constant probability 
density for the sources. After some algebra (see appendix A), one gets: 

dV(r) ^ M ^M-i lr> ^ 

ai ^M sm ( r<J i) i 1 + cos(r<7i)j . (3.4) 

2=1 

Another approach - which can be generalized to non-uniform populations of sources - includes 
Monte Carlo simulations, in which we randomize the coordinates of the candidate sources in 
both RA and dec. Here the Monte Carlo method is used with 10' 5 iterations, averaging the 
distribution of r over the number of iterations, so the resulting histogram is practically free 
from statistical errors associated with the finite number of candidate sources, M. 

When comparing the r-distribution of the data with the null one, the question to be 
answered is how compatible they are, i.e., how likely it is that the former might be a particular 
realization of the latter. To answer this question quantitatively, we use the p-value, defined 
as the probability that the null case produces, in one trial (see sec. 6 for a discussion of 
probabilities for multiple trials), a number of coincidences (r < 1) equal to or larger than 
the one observed in the actual data. Clearly, the larger the excess of the data over the null 
distribution at r < 1, the lower the p- value. Here the p- value is obtained by examining the 
10 5 Monte-Carlo-generated candidate source sets, and finding the percentage of these that 
have a number of coincidences equal to or exceeding the observed one at r < 1. 

4 Here we are using the “nearest neighbor” version of the method, which leads to identifying the nearest 

candidate as the most likely true source of a given neutrino. This may lead to false attributions if more than 
one candidate overlaps with the neutrino data point within the error. Considering the sparseness of the data 
and of the candidate sets we use, we estimate that the chance of false attribution be low. However, for future 
studies with larger samples, one may generalize the current method to include the distance to the second 
closest candidate as an additional statistical variable. 
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Finally, let us comment on the validity of this method and its underlying assumptions: 

(i) no assumptions are made, nor needed, on the spatial distribution of the neutrino data set, 
and on if/how the data correlate spatially with one another. Indeed (see appendix A), the 
main ingredient here is the probability to find a candidate source within the angular error of 
a given neutrino datum. We have verified that our approach is valid for both uniform and 
non-uniform spatial distributions of the data in the sky. 

(ii) We stress that visually comparing the entire r-distribution with the randomized one (as 
shown in figures 1, 2, 3) only has indicative value. This is because the histogram of the 
data is affected by large statistical errors associated with the small number of events in each 
bin and with the small number of candidates in each set. Therefore, we recommend relying 
mostly on the p-value as an indicator of compatibility. 


3.1 Combining neutrino events 

The larger the angular uncertainty of an IceCube neutrino event, the more difficult it becomes 
to disentangle its counterpart and provide a useful statistical evaluation of a population of 
candidates. We note that there are several instances of two or more neutrino shower events 
that are compatible, within the errors, with a common origin in space. It might be of interest, 
then, to explore a common origin hypothesis, and treat the positionally overlapping events 
as different measurements of the position of the assumed same source. In this framework, the 
position of the source can be known more precisely by combining the multiple measurements 
into a single one, using the standard theory of measurement and errors. 

We caution the reader that this hypothesis may imply un-physically large neutrino 
fluxes for individual sources and therefore may be implausible. For this reason, the exercise 
of combining events should be regarded as a useful check, that does not carry the same 
significance as an analysis involving all IceCube-reported neutrino events. We only present 
the results of the “combining” method as a supplement to the more promising findings. 

To combine the neutrino shower events, an iterative procedure is used. At each iter¬ 
ation, all of the N\ possible pairs of data points are considered, with a weighted-distance, 
Cij = Sij/(ai + Uj). If the lowest value of C\j is C \2 < 1 (i.e., the two data points have 
overlapping errors), the corresponding potential pair is combined into a single measurement 
of the position, with its resulting error, as follows: 


OL c — 


Lj=1 a i ~ a i 

v ^2 —2 

Z-^i=l i 



2 



2=1 


(3.5) 


The new neutrino position, (a c ,<5 c ), and error, a c - which is smaller than either of the two 
initial errors, <J\ and 02 - are recorded and replace the original pair of events. 

The process is then repeated, until all overlapping neutrino events have been combined. 
For an example of the “combined” neutrino positions, see the all-sky plots in figures 2 &; 3. 


4 Analysis of possible counterparts 

In this section, the motivation for selected groups of sources and the results obtained are 
presented. Specifically, we give the distribution of the weighted distance, r, (see section 3), the 
p-value, and the excess of events in the first bin (r < 1), AN, relative to the null distribution. 
We provide, as an example, a test consistent with a null distribution to demonstrate a proof 
of method. However, we choose to focus mostly on the most interesting results found for 
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Figure 1. Results for the eleven blazars from the Fermi -LAT 3FGL catalog, for which the 10-100 
GeV gamma-ray flux density is I 7 > 1 x 10~ 9 photons cm -2 s -1 . Left: The equatorial (J2000) 
coordinate sky map of these candidate sources; we distinguish between BL Lac objects (light purple 
diamonds), and flat spectrum radio quasars (FSRQ; green stars). The map also shows the 35 IceCube 
neutrino events as black dots, with their median angular error (pink ellipses). The dot size indicates 
the energy of the neutrino event (see legend). The solid grey line represents the galactic plane. Right: 
The distribution of the weighted-distance to the nearest candidate source of the 35 neutrino events 
(solid, blue). The null distribution, determined via 10 s iterations of a randomized population, is also 
shown (pink with hash marks); purple indicates the overlap of the two histograms. The legend gives 
the excess of the true distribution relative to the null in the the first bin (r > 1), and the p-value. 


groups of starburst galaxies and superbubble regions. A summary of the results for which a 
possible correlation is indicated is given in table 1. 

4.1 “Null” results 

To underline the validity of our method, a test involving the brightest gamma-ray emitting 
blazars is shown in figure 1. A complete description of the selection criteria for the blazar 
candidates can be found in section B. Figure 1 shows, on the left, the location of the eleven 
brightest sources on the sky and on the right, the resulting r distribution of the neutrino 
events. Overall, this distribution is consistent with the null case. Five neutrino data points 
include a blazar within their median error, while six are expected in the null distribution. The 
value p ~0.76 is found for the p- value (see also table 3), leading to the conclusion that there 
is no indication of a causal correlation between the IceCube neutrino events and the brightest 
blazars. This confirms the conclusions found with previous, positionally-blind selections of 
AGN [25, 26], 

A description of the analysis involving additional groups tested for which a “null” result 
was obtained can be found appendix B. Furthermore, a complete summary of the “null” 
results can be found in table 3 of appendix B. 

4.2 Star-forming activity 

Galaxies undergoing star formation at high rates, Rsf — 1CP 1 — 10 2 M 0 yr^ 1 [40], usually 
caused by the disruption or a merger of galaxies, are known as star-forming galaxies. Galaxies 
with star formation rates up to Rsf — 20 M 0 yr -1 , or with typical supernova rates Rsn — 0.3 
yr -1 , and more commonly observed close-by hosting spiral disks, are referred to as starbursts. 
The emission observed over the entire electromagnetic spectrum of star-forming galaxies is 
dominated by the evolutionary processes and environments of stars. A defining feature is 
their luminous infrared emission, peaking just short-ward of 100 /jrn, which is a product of 
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dust absorbing UV radiation from massive stars and re-emitting it in the IR. Star-forming 
galaxies host large populations of objects associated with hadronic gamma-ray emission — 
including supernova remnants (SNR) [41], pulsar wind nebulae (PWN) [42], various types of 
explosions associated with supernova (SN) [43-47], and superbubbles [48] — making them 
prime candidates for CR acceleration and high-energy neutrino emission [9]. Much of this 
activity is connected to the interplay between massive stars and their surrounding media. 

Massive stars (M > 8 Mq) form in dense molecular clouds and live relatively short 
lives (~ 10 6 yrs) before exploding as supernovae. They are typically observed in unbound, 
groups of 0(10— 100) O and B stars, or “OB associations” [49]. The superimposed effects of 
their stellar winds and SN explosions create giant (> 150 pc) cavities of hot, tenuous plasma, 
known as superbubbles [50]. Indeed, about 85% of core-collapse SN occur in superbubbles, 
and starburst galaxies each contain hundreds of these regions [51]. 

Superbubbles are believed to be the origin of “galactic fountains” [52, 53]. In this 
scenario, a star-forming region clears out surrounding gas and dust via stellar winds and SN 
explosions. As this region grows, it will preferentially expand into lower density environments 
and, therefore, in a direction perpendicular to the galactic plane. In time, the star-forming 
region bursts through the galaxy’s disk, exposing the gas and CRs to the halo. The strong 
magnetic fields contribute to this “fountain” effect and direct the collimation of CRs outward. 
This naturally leads to an amplification in CR acceleration [54, 55]. Within the superbubble, 
particle acceleration may be affected by several different processes: shock acceleration in the 
winds, shock acceleration during SN explosions, and second order Fermi processes in the 
turbulent magnetic field deriving from merging stellar winds and SN ejecta [48, 56]. 

Since the detection of high-energy neutrinos, starburst galaxies have been studied as 
promising possible counterparts of IceCube sources [5-8, 17, 27, 43, 44, 57-59]. These dis¬ 
cussions largely rely on the rapid redshift evolution of star-forming galaxies (peaking during 
z ~ 1 — 3 [60]), implying that the majority of a diffuse high-energy neutrino flux should orig¬ 
inate at redshift z > 1. Sources at such large distances are difficult to resolve individually as 
neutrino point-sources, however, due to their low flux and high space density. A positional 
association with two IceCube events has been suggested in a study that identified starburst 
galaxies as the primary candidates to produce the Telescope Array’s UHECR (E > 60 TeV) 
excess [17]. Another interesting connection was the similarity between the gamma-ray flux 
measured by Fermi -LAT within 20° of the Galactic center and the IceCube neutrino flux 
[14]. The proposed scenario suggested increased star-forming activity coupled with CR con¬ 
finement by strong magnetic fields as the most likely origin of the gamma-ray and neutrino 
fluxes. Also within the inner galaxy, the young, compact star cluster Westerlund 1 appears 
to agree well with hard spectra, gamma-ray emission models and with the relative location 
of a number of IceCube neutrino events [61]. 

In the following subsections, we present analyses for groups of (i) starburst galaxies and 
(ii) starburst galaxies plus local star-forming regions and superbubbles. More details on the 
individual candidates are given in section 5. 

4.2.1 Starburst galaxies 

As a selection criterion for starburst galaxies, we made use of their emission at 100 /jrn, which 
is an indicator of star formation (see section 2.2). Specifically, a flux density cut S(100 //m) 
> 250 Jy [32, 33] was imposed, selecting seven starburst galaxies. Two of these, M 82 [62] 
and NGC 253 [63], also have been detected in TeV gamma-ray energies and appear in the 
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Figure 2. Top: The same as figure 1, for all the starburst galaxies that pass the cut on the flux 
at 100 pm, S(100 pin) > 250 Jy [32]. We distinguish those that have been observed by Ferrm-LAT 
(blue squares), and those that have been detected by both Ferrm-LAT and in TeV gamma-rays 
(green triangles). The latter, M 82 and NGC 253, are also the only starburst galaxies detected at 
TeV energies. The remaining starburst galaxies are shown as magenta pentagons. See table 2 for the 
names and coordinates of each starburst. Bottom: Same as the top figures, but for the 25 “combined” 
neutrino events dataset (see section 3.1). The points that are the result of combining two or more 
events are shown as crosses. 


Ferrm-LAT 3FGL catalog (classified as sbg there) as well. Two more of the seven objects 
appear in the Fermi-LAT 3FGL catalog, but have not been seen at TeV energies. 

Figure 2 shows the seven starburst galaxies on the sky map. All of them lie within a 
neutrino event error as reported by IceCube. The r-distribution of the data has two peaks, 
one at 0 < r < 1 and another at 5 < r < 6. This second peak is not expected in the 
null distribution, and might be the result of a statistical fluctuation. The first peak has 
eight neutrino events, an excess AA = 3.8 above the background or null distribution. This 
corresponds to a p-value of 0.042, indicating that there is only a 4% probability that the 
outcome we find is obtained with uniformly distributed sources. This is only moderately 
significant, but sufficiently interesting to motivate further investigation. 

In this spirit, in figure 2 we also show the distribution of our “combined” data set (see 
section 3.1). Five of these data overlap with a starburst galaxy within the positional error, 
constituting an excess of AN = 3.5 relative to null, with a p-value of 0.003. Intuitively, the 
higher significance relative to the uncombined data set can be understood considering that 
the combined positional errors cover a significantly smaller fraction of the sky, so that acci¬ 
dental (i.e., non-causal) coincidence is less likely. The high significance should be contrasted, 
however, with potentially unphysical aspects of the combining exercise. Therefore, the most 
meaningful conclusions here are that, after combining overlapping neutrino data, an excess 
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Figure 3. The same as figure 2, for gamma-ray detected star-forming regions — superbubbles and 
starburst galaxies. In addition to TeV and Fermi- LAT detected starburst galaxies (shown in figure 2), 
this analysis includes the superbubbles detected in TeV emission (down-pointing triangles). One of 
these objects, the Cygnus X Cocoon, was also reported in the 3FGL catalog. 


in the r < 1 persists, while the second peak at r ~ 5 does not, and therefore the former may 
be regarded as more robust. 

In addition to the results shown in figure 2, in table 1 we report further tests on the star 
formation hypothesis. In one test, our method was applied to a reduced set of candidates, 
including only the four 3FGL starburst galaxies. All of these lie within a neutrino event 
error, and the significance of the correlation is in-line with the result for the seven starbursts 
(AN = 3.3 and p = 0.046 for the full data set, and p = 0.001 for the combined one). 

A second test was performed to include a deviation from the hypothesis of uniform 
candidate distribution. Indeed, both gamma-ray and IR observations suffer from source 
confusion and a decreased sensitivity in the direction of the galactic plane [31, 33]. We 
attempted to account for this by excluding all the points that fall within |6| < 10° of the 
galactic plane from the Monte Carlo-generated null distribution. The results of this test 
show a stronger correlation compared to figure 2 for the IceCube neutrino events (p-value of 
0.034), as well as for the “combined” events (p -value of 0.002). In other words, allowing for 
the extragalactic Monte Carlo sources to avoid the galactic plane - as the Fermi and IRAS 
detections do - strengthens the conclusion of figure 2: extragalactic star-forming galaxies are 
a possible neutrino source. 

Additional preliminary tests for which the randomization of the Monte Carlo was per¬ 
formed (i) by varying only the right ascension and not declination and (ii) with a distribution 
equivalent to the sensitivity of IceCube, indicate that the randomizations done in both RA 
and declination produce the most conservative p-values. 
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Candidate 

Catalog(s) 

Selection 

Cand. 

count 

AN 

p -value 



Criteria 

number 

(r < 1) 


(r < 1) 

Starburst 

TeVCat, 3FGL 

starburst 

4 

6 

3.3 

0.046 





[4] 

[3.1] 

[0.001] 

Starburst 

TeVCat, 3FGL 

S(100 pm) > 250 Jy 

7 

8 

3.8 

0.042 


IRAS 100 pm 



[5] 

[3.5] 

[0.003] 

Starburst 

TeVCat, 3FGL 

same as above, 

7 

8 

3.9 

0.034 


IRAS 100 pm 

randomize 
with 6| > 10° 


[5] 

[3.6] 

[0.002] 

Star-form. 

TeVCat, 3FGL 

starburst, 

7 

10 

5.8 

0.003 



superbubble, 
star-form, region 


[6] 

[4.5] 

[<0.001] 


Table 1 . Summary of results presented in figures 2 and 3, and additional tests not shown in graphics. 
For each class of candidates, the table gives the number of objects in the sample (column 4). Entries 
in brackets refer to the combined data set, as explained in section 3.1. Columns 5 and 6 give the 
number of neutrino associations with weighted distance r < 1 from the nearest candidate and the 
excess over the number of associations expected in the null case, AN, respectively. 

4.2.2 Starburst galaxies and local superbubble regions 

Motivated by the interesting result for starburst galaxies of figure 2, let us now focus our 
attention more closely on local star-forming activity. We generated another set of candidate 
sources, with the criterion of supplementing the four gamma-ray detected starburst galaxies 
with superbubbles and star-forming regions that have been detected in TeV gamma rays. This 
addition amounts to two superbubbles and one bubble, or massive star-forming, region.'’ One 
of the superbubbles, the Cygnus Cocoon, is also the only designated star-forming region in 
the 3FGL catalog. 

It was found (figure 3) that ten neutrino events overlap with a candidate within their 
error (r < 1), amounting to an excess of AN = 5.8 and a p -value = 0.003. This suggests a 
low degree of compatibility with the null case. A suggestion of a correlation improves further 
when considering the “combined” neutrino events, where we find 6 positional associations 
with only ~ 1.5 predicted, and the p-value falls below 0.001 (see table 1 for a summary). 

In summary, in section 4 candidate sources were examined related to AGN and/or star 
formation. While the former present a good compatibility with the null case, an interesting 

5 Note that two of the three added candidates are in our own Galaxy (see section 5). Therefore, strictly 
speaking, the set of candidates sources in this subsection can not be a realization of the uniform distribution, 
which is a condition for the validity of the method used here. In this respect, the result obtained with only 
the four starburst galaxies (table 1) can be considered more robust. 
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Name 

RA 

(J2000) 

dec 

(J2000) 

Class 

D l 

[Mpc] 

v ID 

NGC 253 

00 27 34 

-25 17 22 

starburst 

3.1 

7, 10, 21 

NGC 1068 

02 42 43 

-00 01 33 

starburst 

13.7 

1 

[IC 342] 

03 46 49 

±68 05 46 

starburst 

4.6 

31 

30 Dor C 

05 35 55 

-69 11 10 

superbubble 

0.05 

19 

M 82 

09 55 53 

±69 40 46 

starburst 

3.6 

31 

NGC 4945 

13 05 29 

-49 26 03 

starburst 

3.9 

35 

[M 83] 

13 37 01 

-29 51 57 

starburst 

3.6 

16 

W 49 A 

19 10 27 

±09 11 25 

star-form region 

0.011 

25, 33, 34 

Cygnus Cocoon 

20 28 41 

±41 10 12 

superbubble 

0.002 

29, 34 

[NGC 6946] 

20 34 52 

±60 09 13 

starburst 

5.3 

34 


Table 2. The candidate star-forming sources considered in figures 2 and 3, with their equatorial 
coordinates (columns 2 and 3). The names in brackets are those objects that have not been detected 
in gamma-rays, but appear among the brightest starbursts in the IRAS 100 /im catalog (see section 
2.2). Column 5 gives the distances from Earth of each object, taken from [32] for starbursts, [64] for 
30 Dor C, [65] for W 49 A and [66] for the Cygnus Cocoon. Also shown are the ID numbers (from 
[3]) of the neutrino events that have weighted distance r < 1 (eq. (3.3)) for each candidate. 


indication of deviation from it has emerged for star formation. While a robust claim is 
premature, this result is sufficiently interesting to prompt a number of tests to assess the 
plausibility of nearby star-forming sites as being the origin of a subset (~ 3 — 6 data points) 
of the IceCube signal. This is theme of the next section of this paper. 

5 Discussion: star formation as a possible origin of high energy neutrinos 

Let us take a closer look at the seven gamma-ray detected objects rich in star-forming activity 
that have emerged in section 4.2, and apply naturalness considerations to estimate if they 
are plausible contributors to the IceCube neutrino data. 

Table 2 summarizes the main facts of these candidates — and of those that have not 
been detected in gamma rays, which will not be discussed here — and shows the ID number 
of the neutrino data that overlap with these to within the error. Their gamma ray spectra 
are shown in figure 4. Below we describe each object individually. 

M 82 The nearly edge-on, starburst galaxy M 82 (NGC 3034), located ~3.6 Mpc away, 
was the first starburst detected in TeV emission, and perhaps the first direct detection of 
an extragalactic source of hadronic gamma-ray emission [62, 67]. It is the prototypical small 
starburst galaxy with an estimated supernova rate Rsn ~0.1-0.3 yr _1 [68, 69], a gamma- 
ray luminosity L(> GeV) ~ 2 x 10 40 erg s~ 4 [70] and a photon number power-law index of 
T = 2.21 ±0.06 over the 100 MeV - 100 GeV energy bands [31]. Interactions with neighboring 
galaxies, prominently the larger spiral M 81, have spurred star-forming activity, particularly 
in the central regions [71]. 

NGC 253 Thought of as a “twin” of M 82, NGC 253 is similarly located at a distance 
of ~3.1 Mpc with comparable infrared luminosity and spectral distribution and is also seen 
approximately edge-on. The nucleus of NGC 253 contains a very active star-forming region 
150 pc across [72] in which SN occur at a rate Rsn ~ 0-1 yr _1 [73, 74], GeV and TeV 
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gamma-ray emission has been detected, implying a luminosity of L(> GeV) ~ 5.6 x 10 39 erg 
s ” 1 [70], and a power-law fit with a photon index of T = 2.34 ± 0.003 is consistent with no 
spectral break in the gamma-ray emission [75]. X-ray [76-78] and radio [79, 80] observations 
have revealed the presence of a hot diffuse halo resulting from a “disk wind” extending ~9 
kpc from the galactic plane. 

NGC 1068 The most distant gamma-ray detected starburst galaxy, at ~ 13.7 Mpc, NGC 
1068 has the lowest detected 100 /jrn flux, yet its luminosity is more than four times greater 
(Lioop.m ~ 8.6 x 10 24 W Hz -1 ) than the other starburst galaxies selected in our sample 
[81]. Indeed, it is the steep far-IR spectrum and its 100 /xm luminosity that classifies this 
object as a starburst [31, 33]. However, its weak active nucleus, surrounded by a region of 
intense star formation extending ~1 kpc [82], has been widely discussed in the literature as 
the prototypical Seyfert. A study comparing the gamma-ray emission of NGC 1068 detected 
by Fermi -LAT with those of M 82, NGC 253, and NGC 4945 cited that its gamma-ray 
luminosity was too high to be explained only by starburst activity [83]. The best-fit photon 
number power-law index for 100 MeV - 100 GeV energies is T = 2.32 ± 0.10 [31]. 

NGC 4945 Also classified as both a Seyfert II and a starburst galaxy, NGC 4945 is a 
nearly edge-on barred spiral ~3.9 Mpc distant. Unlike NGC 1068, the high energy emission 
detected using Fermi -LAT could be explained only in terms of its starburst activity [83]. 
NGC 4945 is one of the brightest 100 /xm sources with a flux of 1330 Jy [33], only slightly 
fainter than M 82. The best-fit photon number power-law index in the 100 MeV - 100 GeV 
energy bands is T = 2.43 ± 0.07 [31]. A cone-shaped plume extending ~500 pc from the 
nuclear region perpendicular to the disk has been detected in X-ray [84] and Ha [85] and is 
believed to be driven by supernovae. 

Cygnus Cocoon The Cygnus Cocoon is a 50 pc wide star-forming region located in the 
Galaxy 1.5 kpc away [ 66 ] with a combined mass of ~ 8 x 10 6 Mq [ 86 ]. It hosts a collection of 
1500-2000 massive OB stars [87] (e.g Cyg OB2), massive star clusters (e.g. NGC 6910), pul¬ 
sars, SNRs (e.g. 7 Cygni), etc., whose combined effects from stellar winds and SN explosions 
have created a superbubble. Fifteen sources have been detected within the diffuse Cygnus 
Cocoon held in the 3FGL catalog, although the classification of most sources is currently 
unidentified and some may be potentially confused with Galactic diffuse emission [31]. The 
detection of the Cygnus Cocoon, using the Fermi -LAT (1-100 GeV) with a flux L 7 ~ 6 x 10 ” 8 
cnr 2 s ' 1 and the Milagro Gamma-ray Observatory with a flux L 7 ~ 3.5 x 10 " 11 cm ” 2 s ” 1 
centered at ~12 TeV [ 88 , 89], revealed a hard spectrum, most likely of diffuse, interstellar 
origin [ 86 ]. The similarity between the emission morphology and IR/optical observed in¬ 
terstellar structure favors a CR origin, concluding that the Cygnus Cocoon is most likely a 
CR superbubble capable of accelerating CRs up to an estimated 150 TeV [ 86 ]. Additional 
detections at TeV energies have been confirmed by the ARGO-YBJ detector [90]. 

W 49 A Also within the Galaxy, at a distance of 11.4 kpc [65], is the W 49 complex, 
hosting one of the most active and luminous (Ljr > 10' Lq) star-forming regions in the 
Galaxy, W 49 A [91]. W 49 A contains ~30 ultracompact HII regions [92, 93] for a total 
mass of ~ 10 6 Mq [91]. Evidence for multiple expanding shells provided from the radiation 
pressure and/or strong stellar winds of massive stars has been observed (~15 pc scale), as well 
as the remnants of gas ejections on larger scales ( r\_/ 35 x 15 pc 2 ) [94], W 49 A was detected 
in TeV gamma-rays at a 4.4<r significance level, using H.E.S.S and re-analyzed Fermi -LAT 
GeV data [95]. 
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30 Doradus C 30 Dor C is a 100 pc wide superbubble in the Large Megallanic Cloud, 
roughly 50 kpc away [64, 96]. It is luminous in radio, optical [97], X-ray (synchrotron- 
emitting) [98, 99] and TeV gamma-rays [96] with a central temperature of 7.4 x 10 6 K. Al¬ 
though neither leptonic or hadronic originating TeV gamma-ray emission can be ruled out, 
conditions in the superbubble provide evidence for magnetic-field amplification combined 
with turbulence in the hadronic scenario, possibly providing conditions for CRp acceleration 
exceeding energies of 3 PeV [96]. 

Next we will discuss if the physics of these star-forming objects supports strong neutrino 
emission. An immediate test of the promise of an object as a IceCube source is to check if 
it is sufficiently powerful to produce at least one neutrino in IceCube. This can be done 
using the observed gamma-ray spectra (figure 4) under the assumption that neutrino and 
gamma-ray emission escaping the source are similar. This is expected to be true (i) when 
pair production can be neglected, (ii) in most applications of the hadronic model where the 
interstellar medium of the source is transparent to gamma-rays (for Njj < 10 26 cm -2 [100] 
at GeV energies), and (iii) the neutrino flux and spectrum trace the gamma-ray spectrum to 
within a factor rs_/ 2 — 3. 

For a neutrino spectrum of the form E~ 2 the neutrino flux required to produce one event 
in IceCube, is such that E 2 (/)m\ ~ 10~ n ergs cm~ 2 s _1 [19], nearly independently of the 
energy of the specific neutrino event considered, and with a Poissonian error of a factor of a 
few. As a conservative test of the plausibility of a neutrino source candidate, we require that 
its gamma-ray spectrum, power-law-extrapolated to the IceCube energy window (E > 30 
TeV), be within an order of magnitude of <t>m, he. E 2 ^ > E 2 (j) rn i n ~ 10 -12 ergs cm~ 2 s _1 . 
This criterion is similar to that in ref. [19]. 

It appears that M 82 and NGC 253 (top two panels of figure 4) do not pass this 
test, since their spectra decline rapidly with energy and their extrapolations fall below 4> m in 
by at least one order of magnitude. A similar situation is realized for the remaining two 
starburst galaxies, NGC 1068 and NGC 4945. However, here the possibility of a hardening 
of the spectrum above a TeV remains open, due to the absence of observations in this regime. 
Slightly more promising are the spectra for the superbubbles/star-forming regions (bottom of 
figure 4): while 30 Dor C and W 49 A also appear to fail the test, the Cygnus Cocoon passes. 
These results agree overall with past studies of starburst galaxies as high energy neutrino 
sources [70, 101-103], where, although with a strong dependence on the parameters, relatively 
soft neutrino spectra were found, reaching up to 150 TeV [102], Among the galactic objects, 
the Cygnus Cocoon has been predicted to be a detectable neutrino emitter [104-106]. 

We found that the basic energetics test disfavors the simplest scenario of starburst 
galaxies with similar neutrino and gamma-ray spectra. However, it leaves open the possibil¬ 
ity that the gamma-ray flux from starbursts is suppressed compared to the neutrino flux, due 
to unaccounted for absorption affecting the gamma-rays or their parent particles (e.g. ref. 
[107]). The most likely possibility is the case of pair production, where the interaction be¬ 
tween gamma-rays and a radiation field of lower energy photons produces electron-positron 
pairs. This mechanism is responsible for gamma-ray attenuation in connection with the 
extragalactic background light, and it effectively steepens gamma-ray spectra through ab¬ 
sorption and also by redistributing gamma-ray photons to lower energies. The cross section 
for this interaction is maximized when the product of the photon energies is ~ 0.1 MeV 2 . 
For example, for gamma-rays of 100 GeV the interaction is maximized with A ~ 0.1 /im UV 
photons and at 10 TeV energies this interaction is maximized for ~ 10 /im mid-IR photons, 
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Figure 4. Gamma-ray spectra of the candidate star-forming sources (figures 2 and 3, table 2) that 
have been detected in gamma-rays. The solid red line is an estimate of the minimum neutrino flux 
needed to produce one event at IceCube with significant likelihood (see text for details); the width 
of the line represents the energy window probed by IceCube. Top: Fermi -LAT spectra of the four 
starburst classified galaxies in the 3FGL catalog [31], shown as filled symbols. The open symbols 
represent the VERITAS TeV data of M 82 [62] and the H.E.S.S. TeV data of NGC 253 [63]. Bottom: 
Same as top, for the nearby superbubbles and star-forming region. The Cygnus Cocoon TeV spectrum 
is from the ARGO-YBJ detector [90], while for W 49 A, the H.E.S.S and re-analyzed Fermi-L AT 
GeV data are shown [95]. For 30 Dor C, the only data available are H.E.S.S TeV observations [96]. 


both of which are abundant in these galaxies. Recent simulations have shown that gamma- 
rays above E ~ 2 — 5 TeV have opacities above r > 1 in the star-forming galaxy ARP 220, 
M 82, and a similar but smaller effect in NGC 253 [102, 108, 109]. 
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Without relying on gamma ray spectra, it is possible to test for compatibility between 
the subset of n ~ 3 — 4 events that might be due to the objects listed in table 2 — which 
are all local, {d < 15 Mpc) — and the whole set of N = 35 neutrino data at IceCube of 
which B ~ 13 should be due to background (see section 2). The argument is that the local 
contribution to the neutrino flux should have a diffuse counterpart, due to similar objects at 
larger distances that can not be resolved individually. The ratio of local and diffuse fluxes 
can be calculated under the assumption that the local sources are representative of most 
objects in their class, and using basic information on the cosmological evolution of the source 
population. We have estimated this ratio for a non-evolving population of identical objects, 
as well as one evolving like the cosmic star formation rate (see e.g., [110]), and we find that 
up to ~ 2% of the total neutrino flux can be from objects with d < 15 Mpc. In comparison, 
from the data we obtain a ratio / = n/(N — B) ~ 0.14 — 0.18. 

The difference between the two results can be interpreted as disfavoring the hypothesis 
of nearby starburst galaxies as IceCube sources. However, it could also indicate a positive 
fluctuation of the very nearby star formation rate with respect to the cosmological average, 
which has been suggested in connection with measurements of the nearby supernova rate 
[111-114]. Alternatively, it could indicate perhaps a mechanism whereby neutrinos are emit¬ 
ted preferentially along the plane of disk galaxies, leading to an enhanced contribution of 
edge-on galaxies to the total flux. We note in this context that 3 out of the 4 gamma-ray 
detected starburst galaxies in table 2 are seen as edge-on, i.e., with an isophotal axis ratio of 
b/a < 0.35 6 , yet in complete samples of nearby and distant galaxies, on average only ~15% 
of the total galaxies are observed with this orientation [115]. Finally, one should consider 
that exotic properties of the neutrino as a particle could increase the ratio /; for example, 
such an effect is predicted in models where a new neutrino interaction with a light mass 
mediator causes absorption of the cosmological contribution to the neutrino flux [116]. 

6 Summary and conclusion 

We addressed the question of a possible contribution to the IceCube data from local objects 
with powerful high energy emissions due to AGN and star-forming activities. Specifically, we 
tested for a statistical positional correlation between the neutrino data and candidate sources 
from different catalogs, selecting among objects of the same class/morphology by imposing a 
minimum flux in either gamma-rays or 100 fim emission. The statistical quantity of interest is 
the weighted distance, r, between a neutrino data point and its nearest candidate source. We 
tested the distribution of r against the “null” hypothesis of random, non-causal, positional 
overlap with uniformly selected candidates sources. 

Results are consistent with a null hypothesis for control samples of nearby stars and 
exoplanets, as well as for blazars and Seyfert galaxies brightest in gamma-ray emission as 
seen in the Fermi- LAT 3FGL catalog (see figure 1 and appendix B). 

More interesting results are found for objects with high star formation rates, such as star- 
burst galaxies, superbubbles and massive star-forming regions. The most significant excess 
at r < 1 is found for the set of seven gamma ray-detected (from TeVCat and 3FGL catalogs) 
star-forming objects, including four starburst galaxies, two superbubbles (one galactic, the 
Cygnus Cocoon) and one galactic star-forming region. Ten neutrino events overlap with a 
candidate within the error. The probability of this occurring in the null hypothesis for one 
trial (p-value) is p = 0.003, as shown in figure 3. 

6 NASA/IPAC Extragalactic Database; http://ned.ipac.caltech.edu/ 
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Similar excesses, although less significant, are seen for different selection criteria, which 
correspond to sets of candidate sources that partially overlap with one another. In particular, 
the p-value is at the level of a few percent for the seven starburst galaxies that are brightest 
in far-IR 100 /im emission (see figure 2). This set includes the same four starburst galaxies 
as in the previous case. If only the latter are considered, the p-value remains consistent with 
p ~ 0.05. 

Taking a more global look at our results, one may wonder about the combined proba¬ 
bility to obtain, in the null case, all the results shown here. Statistically, the K = 8 cases 
analyzed here (tables 1 and 3) represent trials, each of them with a p- value, pi (i = 1,2,..., K). 
Note that only J = 4 trials (stars, blazars, Seyfert galaxies and starburst galaxies) are truly 
independent. Let us focus on the minimum p- value obtained, p m in = 0.003. A commonly 
used quantity [117] is the probability to obtain at least one outcome with p-value p rnin -, over 
the number of trials done. Restricted to the number of independent trials, one would get 
[118] P = 1 — (1— Pmin) J — Jpmin = 0.012. P is frequently referred to as the “post-trial” 
p-value. For the case at hand, of K non-independent trials, a correct estimation of P is 
more complicated and is beyond the scope of our paper. However, generalizing the reasoning 
above, we can be confident that Jp m i n ^ P ^ Kp m i n , i.e., 0.012 P <, 0.024. 

At this point, considerations on the significance of the excess of positional associations 
are necessarily mixed. On one hand, the excess is robust, since it appears for different 
candidate selection criteria and different models of the null case (uniform distribution or 
uniform with subtracted Galactic plane, see table 1). In contrast, however, this result is not 
fully supported by a basic energetic test of neutrino emission tracing gamma-ray emission 
and the assumed redshift evolution of the star-forming population. 

For the seven gamma-ray detected candidates, the gamma-ray spectra were examined. 
Under the assumption that neutrino spectra trace the gamma-ray ones, only one object, the 
Cygnus Cocoon in our galaxy, was found to be sufficiently powerful to produce one event in 
IceCube with substantial likelihood. The remaining candidates can be reconciled with the 
observed excess in the hypothesis that their gamma ray flux might be substantially dimmer 
than the neutrino flux, due to e.g., high absorption. It is also worth noting that the excess 
corresponds to a local contribution (from sources within ~15 Mpc distance) to the total 
neutrino flux of about 10-20%. This is higher than the expected few percent when compared 
to the cosmological evolution of star formation, and could be explained by a local fluctuation 
in star-formation rate relative to the cosmic average or possibly, a preferential orientation for 
neutrino production in disk galaxies. 

Looking ahead, we expect that the situation will become clearer with better statistics at 
IceCube. An update with 54 data points is upcoming (see [119] for preliminary presentations 
at conferences), and subsequent updates are expected at a rate of roughly 12 new events per 
year. Therefore it is likely that, within a year or two, the excess we observe might either 
become disfavored, or confirmed with higher significance. A positive result would have the 
character of discovery, and would be very fertile of theoretical developments on the physics of 
starburst galaxies and other star-forming regions. Questions to be investigated will be how 
opaque these star-forming regions could be to gamma-rays, and if the observational viewing 
angle of star-forming galaxies can affect the measurable neutrino flux. 
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A Formalism: the null hypothesis 


Here we derive the distribution of the weighted distance, r, for the null hypothesis, eq. (3.4). 
Let us consider a population of M candidate sources, uniformly distributed in the sky, so 
that the probability to find a candidate in a unit of solid angle is dp/dQ = 1/(47t). Consider 
now a generic point in the sky, with 6 the angular distance from it. The probability to find 
a candidate at angular distance between 9 and 9 + d9 is, then: 

dp(0) = ^ sin 6d0 . (A.l) 

By integration, one gets the probability to find a candidate at distance larger than 9: 

q(°) = 7,(1 +cos9) . (A.2) 

From these, we can obtain the probability that the nearest source is at angular distance 
between 9 and 9 + d9. This is given by the probability that one source is between 9 and 
9 + d9, and all the other M — 1 candidates are at a larger distance [38]: 

dP(9) = sin 0(1 + cos 9) M ~ 1 d9 , (A.3) 

where the factor M in the numerator is found by the assumption of identical sources. 

Some observations on the distribution dP/d9 in eq. (A.3): 


• dP/d9 = 0 at 9 = 0 and at 9 = n, as expected. It has a maximum at 9 = 9 max , with 


cos 9 max — 1 


(A.4) 


which agrees with the intuition that, for larger M, the most likely distance to the 
nearest source is smaller. In the approximation 9 <C 1, eq. (A.4) gives 9max ~ y 7 2/M. 
The dependence on M -1 / 2 is expected, considering that for M objects populating a 
two-dimensional space the area occupied by each object scales like 1/M. 


• The average of the distribution dP/d9 is found to be, in the limit M 3> 1: 


( 9 )^ 



(A.5) 


It is useful to express the distribution in eq. (A.3) in terms of a weighted distance, 
r = 9/a, with a a constant (in our specific application, a is the angular error on the measured 
neutrino position, see section 3): 

dP(r, a) M M _ x 

— 7Z — = sm ( rcr ) l 1 + cos(rcr)J . (A.6) 
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For a set of N neutrino data, each with error ai (i = 1,2,., iV), the distribution of r is the 

sum: 

dp ( r ) _ f dp {r, ai) , A 

dr ^ dr 1 ' ’ 

i= 1 

The combination of eqs. (A.7) and (A. 6) gives the expression in eq. (3.4). We have checked 
that this result coincides with the Monte Carlo-simulated one for a uniformly distributed 
population of candidates. 

B Null results 

In this section, we present a statistical analysis of candidate sources tested that have produced 
results consistent with a null distribution. We choose two types of sources, the brightest stars 
and closest exoplanets, expected not to correlate with the location of IceCube neutrino events, 
to show a confirmation of the validity of our method. We additionally present the results 
of the brightest emitting gamma-ray blazars and Seyfert galaxies. As done in section 4, we 
present the distribution of the weighted distance, r (see section 3), the p- value, and the excess 
of events in the first bin (r < 1), AN, relative to the null distribution. A summary of null 
results is in table 3. 

Stars We selected the seven brightest (in my) stars as seen in Smithsonian Astrophysical 
Observatory (SAO) Star Catalog [120], therefore having an apparent V magnitude of V < 0.3. 
Seven were chosen to ease computation and compare directly with the star-forming sources. 
It was found that only two neutrino events overlapped with a stellar source resulting in a 
high likelihood, a p-value of 0.977, that the positional coincidences of the brightest stars are 
consistent with the null distribution (see table 3). 

Exoplanets As an additional control sample, we found the seven closest stellar systems 
hosting exoplanets. To compile the “candidate” sources, confirmed exoplanet systems from 
among the Exoplanet Orbit Database (EOD; [121]) were chosen to have a host star distance 
of d < 4.95 pc. Only one candidate source was considered for systems with more than 
one confirmed exoplanet. In this case, five neutrino events overlapped with a stellar system 
hosting an exoplanet, showing an excess of AN = 0.3 over the expected 4.3 associations in 
the null distribution at r < 1, resulting in a p = 0.530 (see table 3). 

Blazars (AGN) Blazars are types of actively accreting AGN whose variable emission 
largely dominates their hosting galaxy, and for which highly relativistic beams are oriented 
along the line-of-sight [122, 123]. They are divided into two major classes. Those displaying 
strong and broad optical emission lines are usually called flat-spectrum radio quasars (FS- 
RQs), while objects with no broad emission lines are called BL Lacertae (BL Lac) objects 

[124] . Blazars are considered to be natural mechanisms for high-energy particle acceleration. 
In particular, the acceleration of protons may explain: (i) the energy transfer from the central 
engine over distances as large as 1 pc, (ii) the heating of a dusty disc in the nucleus over dis¬ 
tances of several 100 pc, and (iii) a near-infrared cut-off of the synchrotron emission in the jet 

[125] . However, recent literature [34, 35] suggests that the currently favored mechanism for 
driving the high energy emission from blazars is a population of electrons accelerated to TeV 
energies, typically through Fermi acceleration by shocks in the AGN jet. TeV gamnra-ray 
emission results from inverse Compton scattering off relativistic electrons, and the electrons 
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Candidate 

Catalog(s) 

Selection 

Criteria 

Cand. 

number 

count 
(r < 1) 

AN 

p -value 
(r < 1) 

Blazar 

3FGL 

-^10—100GeV > 

11 

5 

-1.0 

0.764 



10 -9 ph. cm -2 s _1 


[1] 

[-1.2] 

[0.938] 

Seyfert 

3FGL 

Seyfert I & II 

6 

6 

2.2 

0.165 





[2] 

[0.7] 

[0.368] 

Stars 

SAO Star Cat 

V < 0.3 

7 

2 

-2.7 

0.977 





[1] 

[-0.5] 

[0.828] 

Exoplanets 

EOD 

d < 4.95 pc 

7 

5 

0.3 

0.530 





[1] 

[-0.5] 

[0.828] 


Table 3. The same as in table 1 but for the “null” results presented in figures 1, 5, and from 
appendix B not shown in graphics. 


cool by radiating X-ray synchrotron emission. The strong correlation often observed between 
X-ray and TeV gamma-rays from blazars indicates a possibly common origin. 

In our analysis, we searched for the brightest AGN-classified objects in the 3FGL cat¬ 
alog. This included sources such as blazars (both BL Lac &; FSRQ), AGN, radio galaxies, 
unidentified blazar candidates, and quasars — or those classified as BLL, bll, FSRQ, fsrq, 
agn, RDG, rdg, BCU, bcu, ssrq in 3FGL. A lower bound on the photon flux was imposed: 
L 7 > 10 -9 photons cnh’ s -1 in the 10-100 GeV band. Out of the 11 objects passing this 
criterium, 10 are BL Lac blazars and 1 is of FSRQ type. Indeed, a higher flux in gamma-ray 
emission is expected when the orientation of the AGN jet is pointed in the direction of our 
viewing angle (blazars) compared with larger jet orientation angles (radio galaxies, AGN 
classified, quasars, etc.) [34], Figure 1 shows the resulting r distribution of the neutrino 
events. Overall, this distribution is consistent with the null case. Five neutrino data points 
include a blazar within their median error, while six are expected in the null distribution. The 
value p ~0.76 is found for the p-value, leading to the conclusion that there is no indication 
of a causal correlation between the IceCube neutrino events and the brightest blazars. This 
confirms the conclusions found with previous, positionally-blind selections of AGN [25, 26]. 

Seyfert galaxies Seyfert classified galaxies are characterized by a bright nucleus, with an 
AGN strength in emission that is below that of a quasar or blazar L < 10 11 Lq 7 [126, 
127]. Seyfert galaxies typically have spiral morphologies and active areas of star formation 
surrounding the nucleus. Their emission has contributions from both the central AGN and 
star-forming activity in the galaxy disk. In the picture of AGN unification [123], Seyferts 
are identified as Type II when viewed edge on and Type I when the jet is oriented along the 

' Here Lq is the luminosity of the sun, Lq = 3.846 x 10 33 erg s 1 . 
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Figure 5. The same graphics as figure 1, now for the six Seyfert classified objects detected in the 
3FGL catalog, divided in Seyfert I objects (magenta left-pointing triangles) and Seyfert II objects 
(blue right-pointing triangles). 


line-of-sight. Moreover, these AGN are relatively abundant at low z, as they are thought to 
be the evolutionary by-products of quieting quasars and blazars at higher redshifts [128]. 

For our analysis, we chose all six Seyfert classified objects in the Fermi- LAT 3FGL 
catalog. Of these, five are Seyfert I’s (designated NLSY1 or nlsyl in the 3FGL catalog) and 
one is Seyfert II (classified as sty in the 3FGL catalog). Results from this analysis are shown 
in figure 5. Three neutrinos events overlap with a candidate within the error. The excess in 
the first bin of the r-distribution is AN ~ 2.2, and a p -value of 0.165 is found. Hence, Seyferts 
do not constitute a significant signal; however, our results suggest that Seyfert galaxies may 
warrant further investigation. 
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